Average Bubble Profile

import numpy as np
import math
import matplotlib
#matplotlib.rcParams.update({'font.size': 12})
from matplotlib import gridspec
import matplotlib.pyplot as plt
from collections import deque
import scipy as scp
import scipy.optimize as sco
import scipy.signal as scs
import scipy.special as ssp
import scipy.integrate as sci
import scipy.interpolate as intp
import statistics as stat
import random
from functools import partial
from operator import eq
from itertools import zip_longest, compress, count, islice, groupby, cycle
from labellines import labelLine, labelLines
import os

np_load_old = np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)

Simulation Parameters

nLat = 4096
nSims = 100
minSim = 0
simstep = 1

nu = 2.*10**(-3)
lamb = 1.5; print('lamb = ', lamb)
m2eff = 4. * nu * (- 1. + lamb**2); print('m2eff = ', m2eff)
lenLat = 2 * 50. / np.sqrt(2. * nu); print('lenLat = ', lenLat)
phi0 = 2. * np.pi / 4.; print('phi0 = ', phi0)

alpha = 8.
nCols = 1
phi_initial = np.pi
dt_phi_initial = 0.
mask = 4*phi_initial

############################################################
nyq = nLat//2+1
spec = nyq
hLat = nLat//2
dx = lenLat/nLat
dk = 2.*np.pi/lenLat
dt = dx/alpha
dtout = dt*alpha
light_cone = dtout/dx

lamb =  1.5
m2eff =  0.01
lenLat =  1581.1388300841897
phi0 =  1.5707963267948966

titles = [r'$\phi(x)$', r'$\partial_t \phi(x)$', r'$|\nabla \phi(x)|^2$', r'$V(\phi(x))$']
plots_file = '/home/dpirvu/big_plot_file/thick_wall_average_bubble/'
pickle_location = '/home/dpirvu/pickle_location/thick_wall_average_bubble/'
cutout_location = '/gpfs/dpirvu/thick_wall_average_bubble/'
suffix = '_for_phi0{:.4f}'.format(phi0)+'_lambda'+str(lamb)+'_len{:.4f}'.format(lenLat)+'_x'+str(nLat)

def bubbles_file(min, max):
    return cutout_location+'bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def centred_bubbles_file(min, max):
    return pickle_location+'centred_bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def bubbles_at_rest_file(sim):
    return cutout_location+'restbubble_sim'+str(sim)+suffix+'.npy'
def average_bubble(min, max):
    return pickle_location+'average_bubble_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def average_of_N_bubbles(Nbubbles):
    return pickle_location+'average_of_'+str(Nbubbles)+'_bubbles'+suffix+'.npy'

instanton_location = '/home/dpirvu/inst/instantons/dev/thick_wall_instanton_sim.dat' #/np.sqrt(m2eff)/dx
a = np.genfromtxt(instanton_location)
coleman_profile = np.pi+a[:,0]
xoffset = nLat//8
temp = coleman_profile[len(coleman_profile)//2-xoffset:len(coleman_profile)//2+xoffset]
peaks, _ = scs.find_peaks(temp)
fwhm, height, left_ips, right_ips = scs.peak_widths(temp, peaks, rel_height=0.5)
radius_Coleman_bubble = max(fwhm)/2.
filter_size = radius_Coleman_bubble*dx
print('filter_size = ', filter_size)

plt.plot(np.arange(len(temp))-xoffset, temp)
plt.plot(peaks-xoffset, temp[peaks], "x")
[plt.hlines(height[i], left_ips[i]-xoffset, right_ips[i]-xoffset, color="C2") for i in range(len(fwhm))]
plt.xlabel(r'$r_{\rm E}$'); plt.ylabel(r'$\phi$')

filter_size =  19.550844379806723

Text(0, 0.5, '$\\phi$')

def V(phi):
    return ( -np.cos(phi) + 0.5 * lamb**2. * np.sin(phi)**2. ) * 4. * nu
def dV(phi):
    return ( np.sin(phi) + 0.5 * lamb**2. * np.sin(2.*phi) ) * 4. * nu

right_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[np.pi, 2*np.pi], method='bounded')
left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[0, np.pi], method='bounded')
right_left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[2*np.pi, 3*np.pi], method='bounded')

def F(x):
    return V(x) - V(phi_initial)
phi_upper_bound = sco.fsolve(F, 5)[0]
phi_lower_bound = sco.fsolve(F, 1)[0]
phi_upper_lower_bound = sco.fsolve(F, 8)[0]

############################################################
wall_tension, err = sci.quad(lambda x: np.sqrt(2*(V(x) - V(2*np.pi))), np.pi, 2*np.pi); print('wall_tension = ', wall_tension)
epsilon = V(np.pi) - V(2*np.pi); print('epsilon = ', epsilon)
R_coleman = wall_tension/epsilon; print('R_coleman = ', R_coleman)
wall_thkn = right_phi_at_V_max.x / np.sqrt( V(right_phi_at_V_max.x) - V(np.pi) ); print('wall_thkn = ', wall_thkn)
############################################################

matplotlib.rcParams.update({'font.size': 12})
fig = plt.figure(figsize=(7, 4))
plt.plot([i for i in np.arange(-np.pi/2, 5*np.pi/2, np.pi/100)], [V(i) for i in np.arange(-np.pi/2, 5*np.pi/2, np.pi/100)])
plt.plot(phi_upper_bound, V(phi_upper_bound), 'ro')
plt.plot(phi_lower_bound, V(phi_lower_bound), 'co')
plt.plot(phi_upper_lower_bound, V(phi_upper_lower_bound), 'yo')
plt.xlabel(r'$\phi(x)$'); plt.ylabel(r'$V(\phi(x))$'); plt.show()

wall_tension =  0.46497413771169294
epsilon =  0.016
R_coleman =  29.06088360698081
wall_thkn =  80.67292631362025

VISUALISING SIMULATION DATA

all_data, sims_to_keep = np.load(bubbles_file(0, nSims))
normal = [phi_initial, np.mean([sim[1][0] for sim in all_data]), np.mean([sim[2][0] for sim in all_data]), V(phi_initial), 0.5*np.mean([sim[1][0] for sim in all_data])**2]
#normal = [phi_initial]

def add_mask(field_slice, threshold):
    return field_slice * [0 if np.abs(i) > threshold else 1 for i in field_slice]

def smoothen(field_slice, sigma, bool):
    N = len(field_slice)
    if bool:
        pbc = [j-1 if j < N/2+1 else N-j+1 for j in range(1, N+1)]
    else:
        pbc = [j-1 for j in range(1, N+1)]
    window = [np.exp(- 0.5 * (x*dx/sigma)**2) / np.sqrt(2*np.pi) / (sigma/dx) for x in pbc]
    window = window / sum(window)
    spectral_filter = np.fft.fft(window, len(field_slice))

    fft_field_slice = np.fft.fft(field_slice, len(field_slice))
    smooth_fft_field_slice = [spectral_filter[k] * fft_field_slice[k] for k in range(len(fft_field_slice))]
    smooth_field_slice = np.fft.ifft(smooth_fft_field_slice, len(field_slice))
    return np.asarray([i.real for i in smooth_field_slice])

def find_slice_peaks(field_slice, peak_threshold):
    """ Finds x coordinate of peaks in masked field with mask applied at threshold. """
    peak_coord = signal.find_peaks(field_slice, height = peak_threshold)[0].tolist()
    if field_slice[-1] > peak_threshold and field_slice[0] > peak_threshold and field_slice[-1] != field_slice[0]:
        if field_slice[0] > field_slice[-1] and field_slice[0] > field_slice[1]:
            peak_coord.append(0)
        elif field_slice[0] < field_slice[-1] and field_slice[-1] > field_slice[-2]:
            peak_coord.append(len(field_slice)-1) # this minds potential boundary discontinuities
    peak_heights = [field_slice[coord] for coord in peak_coord]
    return peak_coord, field_slice.tolist().index(np.max(peak_heights))

def truncateNum(num, decimal_places):
    StrNum = str(num)
    p = StrNum.find(".") + 1 + decimal_places
    return float(StrNum[0:p])

def time_at_fraction(bubble, frac, limit):
    T, N = len(bubble), len(bubble[0])
    right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
    time_list = [t if (right_phi_x[t] <= N*frac) else 0 for t in range(T)]
    return next((t for t in time_list[::-1] if t != 0), 0)

def time_at_size(bubble, size, limit):
    T, N = len(bubble), len(bubble[0])
    right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
    time_list = [t if (right_phi_x[t] <= size) else 0 for t in range(T)]
    return next((t for t in time_list[::-1] if t != 0), T-1)

def center_bubble(bubble):
    bubble0 = bubble[0]
    #truncate bubble where it expands relativistically
    limit = phi_upper_bound
    tdecap = time_at_fraction(bubble0, 0.6, limit)
    T, N = len(bubble0), len(bubble0[0])

    tcheck = time_at_fraction(bubble0, 0.02, 2*phi_initial)
    slice = smoothen(bubble0[tcheck], dx, True)
    xmax = slice.tolist().index(max(slice))
    angle = int(N//2) - xmax # rotation angle needed
    print(angle)
    bubble = np.asarray([[np.roll(bubble[col][t], angle) for t in range(tdecap)] for col in range(len(bubble))])

    # check which sides the COM travels
    tcheck = time_at_fraction(bubble[0], 0.02, 2*phi_initial)
    slice = smoothen(bubble[0][tcheck], dx, True)
    xmax = slice.tolist().index(max(slice))
    xlist = np.arange(nLat//2-500, xmax+500)%nLat
    fld = [[bubble[0][t][x] for x in xlist] for t in range(tcheck-100, tcheck, 2)]
    tv = [np.nanmean([x for x in range(len(fld[0])) if fld[t][x] > np.ceil(limit)]) for t in range(len(fld))]
    decr = np.sum([1 for i, j in zip(tv, tv[1::]) if i >= j])
    incr = np.sum([1 for i, j in zip(tv, tv[1::]) if i < j])
    print('tcheckm, xmax, incr, decr ', tcheck, xmax, incr, decr)
    return bubble, decr - incr #np.nanmean(tv) - N//2

def multiply_bubble(bubble, dir, fold):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    bubble = [np.tile(bubble[col], fold) for col in range(len(bubble))] # multiplies bubbles so tail is kept without pbc
    TT, NN = len(bubble[0]), len(bubble[0][0])
    for t in range(TT):
        a, b = int((TT-t)*light_cone) + N, int((TT-t)*light_cone/2) - N//4
        x1, x2 = np.arange(a, NN), np.arange(b)
        if dir < 0:
            x1, x2 = x1 - a, x2 - (b-NN)
        for x in np.append(x1, x2):
            if 0 <= x < NN:
                bubble[0][t][x] = phi_initial
    return np.asarray(bubble)

#bubble_data = []
for sim in range(len(all_data))[:1:]:
    print(sim)
    col = 0
    simulation = all_data[sim]
    centered, dir = center_bubble(simulation)
    mult = multiply_bubble(centered, dir, 2)
    fig, ax = plt.subplots(1, 3, figsize = (12, 3))
    im0 = ax[0].imshow(simulation[col], aspect='auto', interpolation='none', origin='lower')
    im1 = ax[1].imshow(centered[col], aspect='auto', interpolation='none', origin='lower')
    im2 = ax[2].imshow(mult[col], aspect='auto', interpolation='none', origin='lower')
    clb = plt.colorbar(im0, ax = ax[0]); clb = plt.colorbar(im1, ax = ax[1]); clb = plt.colorbar(im2, ax = ax[2]); plt.show()
#    bubble_data.append(mult)
#np.save(centred_bubbles_file(minSim, nSims), bubble_data)
#check_multiplication = np.load(centred_bubbles_file(minSim, nSims))
#for pic in check_multiplication:
#    plt.figure(0); plt.imshow(pic, aspect='auto', interpolation='none', origin='lower'); plt.show()

0
1385
tcheckm, xmax, incr, decr  914 2048 20 29

def rapidity(vel):
    return np.arctanh(vel)
def gamma(v):
    return ( 1. - v**2. )**(-0.5)

def get_bubble_radius(bubble_slice, filter):
    bubble_slice = smoothen(bubble_slice, filter, True)
    results_half = signal.peak_widths(bubble_slice, signal.find_peaks(bubble_slice)[0], rel_height=0.5)
    return max(results_half[0])/2.

def tanh_pos(x, r0L, r0R, dr, vL, vR):
    wL, wR = dr/gamma(vL), dr/gamma(vR)
    return ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi
def tanh_neg(x, r0L, r0R, dr, vL, vR):
    wL, wR = dr/gamma(vL), dr/gamma(vR)
    return - ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi

def tanh_fit(bubble_slice, axis, prior):
#    plt.plot(axis, bubble_slice, 'ro', axis, [tanh(r, *best_tanh) for r in axis], 'g'); plt.show()
    bounds = ((axis[0], 0, 0, 0, 0), (0, axis[-1], axis[-1], 1, 1))
    if prior is not None:
        if np.mean(bubble_slice) > phi_initial:
            return sco.curve_fit(tanh_pos, axis, bubble_slice, p0=prior, bounds=bounds)[0]
        else:
            return sco.curve_fit(tanh_neg, axis, bubble_slice, p0=prior, bounds=bounds)[0]            
    else:
        if np.mean(bubble_slice) > phi_initial:
            return sco.curve_fit(tanh_pos, axis, bubble_slice, bounds=bounds)[0]
        else:
            return sco.curve_fit(tanh_neg, axis, bubble_slice, bounds=bounds)[0]

def hyperbola1(t, a, b, c):
    return np.sqrt(c + (t - b)**2) + a
def hyperbola2(t, d, e, f):
    return - np.sqrt(f + (t - e)**2) + d

def hypfit_test(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
    try:
        if rr[0] <= rr[-1]:
            fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
            fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
        else:
            fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
            fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
        return True
    except (RuntimeError, TypeError):
        return False
    
def hypfit(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
    if rr[0] <= rr[-1]:
        fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
        fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
        traj1, traj2 = hyperbola1(tt, *fit1), hyperbola2(tt, *fit2)
    else:
        fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
        fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
        traj1, traj2 = hyperbola2(tt, *fit1), hyperbola1(tt, *fit2)
    print('hyperbolic trajectories fit: ', fit1, fit2)
    plt.plot(rr, tt, 'g-', qq, tt, 'r-', traj1, tt, 'y:', traj2, tt, 'b:') # plot the equation using the fitted parameters
    return tt, traj1, traj2

def get_velocities(tt, rrwallfit, llwallfit):
    #uu = wall travelling along with the com; vv = wall travelling against the com; aa = centre of mass velocity; bb = instantaneous velocity of wall
    uu = intp.splev(tt, intp.splrep(tt, rrwallfit), der=1)
    vv = intp.splev(tt, intp.splrep(tt, llwallfit), der=1)

    for i in range(len(uu)):
        if np.abs(uu[i]) > 1:
            uu[i] = np.sign(uu[i])*(1-np.abs(np.abs(uu[i])-1))
        if np.abs(vv[i]) > 1:
            vv[i] = np.sign(vv[i])*(1-np.abs(np.abs(vv[i])-1))
    for i in range(len(uu)):
        if str(uu[i]) == 'nan' :
            uu[i] = - np.sign(vv[-1])*0.999
        if str(vv[i]) == 'nan' :
            vv[i] = - np.sign(uu[-1])*0.999

    aa = ( 1 + uu*vv - np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / ( uu + vv)
    bb = (-1 + uu*vv + np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / (-uu + vv)
    return uu, vv, aa, bb#, ch1, ch2

def velocity(bubble, bool):
    T, N = len(bubble), len(bubble[0])
    limit = phi_upper_bound #phi_upper_bound# - 0.5*np.abs(phi_upper_bound-np.floor(phi_upper_bound))
    prior = None
    data_list = []
    err = 0

    # find where fit can begin
    window = int(N//10)
    if window > 100: window = 100

    tf = T-1
    endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
    if len(endSlice) == 0:
        return 'next'
    while endSlice[0] - window < 0 or endSlice[-1] + window >= N:
        tf -= 2
        endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
        if len(endSlice) == 0:
            return 'next'

    for t in range(tf, 0, -1):
        if len(data_list) == 0:
            peaks = endSlice
            target = int(np.round(np.nanmean(peaks)))
            coord_list = np.arange(peaks[0] - window, peaks[-1] + window)
        else:
            xrange = np.arange(data_list[-1][0]-window, data_list[-1][0]+window+1)
            if all(0 <= i < N for i in xrange) and any(np.cos(bubble[t][i]) > np.cos(limit) for i in xrange):
                peaks = [i for i in xrange if np.cos(bubble[t][i]) > np.cos(limit)]
            else:
                break
            target = int(np.round(np.nanmean(peaks)))
            coord_list = np.arange(target - int(np.round(np.abs(data_list[-1][1]))) - window, target + int(np.round(np.abs(data_list[-1][2]))) + window)

        coords = np.asarray([(bubble[t][i%N], int(i-target)) for i in coord_list])
        if err < 5:
            try:
                r0L, r0R, dr, vL, vR = tanh_fit(coords[:,0], coords[:,1], prior)
                prior = r0L, r0R, dr, vL, vR
                data_list.append([target, r0L, r0R, int(t)])
            except (RuntimeError, ValueError, TypeError):
                prior = None
                data_list.append([target, r0L + np.sign(data_list[-1][1]-data_list[0][1])*light_cone*(T/N), r0R + np.sign(data_list[-1][2]-data_list[0][2])*light_cone*(T/N), int(t)])
                err += 1
        else:
            break
    data_list = np.asarray(data_list[::-1])
    targets, r0Ls, r0Rs, time_list = data_list[:,0], data_list[:,1], data_list[:,2], data_list[:,-1]
    print('tmin, tmax', time_list[0], time_list[-1])

    # get direction of bubble
    radius_diff = np.mean([np.abs(data_list[i,1]) - np.abs(data_list[i,2]) for i in range(len(data_list)//4)])
    if radius_diff > 0:  # if average difference in radius is positive then left radius is on average higher than right radius i.e. bubble travels to the right
        rr, ll = targets + r0Rs, targets + r0Ls
    else:
        rr, ll = targets + r0Ls, targets + r0Rs

    if bool:
        fig, ax0 = plt.subplots(1, 1, figsize = (5, 4))
        ax0.plot(rr, time_list, 'b', ll, time_list, 'y', linewidth='3')
        ax0.set_xlabel('x'); ax0.set_ylabel('t')

    # get velocities from derivative of best fit to wall trajectory
    trunc = 0
    # save copies
    time_list_copy, rr_copy, ll_copy = time_list, rr, ll
    # as lons as it needs..
    while(True):
        # try to fit walls
        if hypfit_test(time_list_copy, rr_copy, ll_copy):
            time_list, rrwallfit, llwallfit = hypfit(time_list_copy, rr_copy, ll_copy)
            uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)
            break
        # if fit fails, shorten walls
        else:
            trunc += 25
            if len(time_list_copy) > 100 + trunc*2:
                time_list_copy, rr_copy, ll_copy = time_list_copy[trunc:-trunc:], rr_copy[trunc:-trunc:], ll_copy[trunc:-trunc:]
            if trunc > 200 or len(time_list_copy) < 100:
                return 'next'
            continue

#    if bool:
#        fig, [ax0, ax1] = plt.subplots(1, 2, figsize = (15, 4))
#        ax0.plot(rr[-len(time_list):], time_list, 'b', ll[-len(time_list):], time_list, 'y', rrwallfit, time_list, 'r:', llwallfit, time_list, 'r:', linewidth='3')
#        ax0.set_xlabel('x'); ax0.set_ylabel('t')
#        ax1.plot(time_list, uu, 'b:', time_list, vv, 'y:')
#        ax1.plot(time_list, aa, 'r', label=r'v COM')
#        ax1.plot(time_list, bb, 'g', label=r'v walls')
#        ax1.axhline(0, color='darkgray', ls=':')
#        ax1.set_xlabel('t'); ax1.set_ylabel('v(t)'); ax1.legend(); plt.show()
#        plt.savefig(plots_file+'velocity_profile'+str(random.randrange(100))+suffix+'.png');

    list = np.abs(aa - bb)
    if any(str(i) != 'nan' for i in list):
        return -aa[list.tolist().index(np.nanmin(list))]
    else:
        return 'next'

def fold(beta):
    return 2

def add_velocities(v1,v2):
    return (v1 + v2) / (1. + v1*v2)

def cut_out_bubble(bubble):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    ext = int(N//8)
    if ext < 100: ext = 200
    xmin = next((i for i in range(N) if bubble0[-1][i] > right_phi_at_V_max.x), 0) - ext
    xmax = next((i for i in range(N)[::-1] if bubble0[-1][i] > right_phi_at_V_max.x), N-1) + ext
    if xmin < 0: xmin = 0
    if xmax >= N: xmax = N-1
    t, tmin = 0, 0
    while not any(bubble0[t][i] > right_phi_at_V_max.x for i in range(xmin, xmax+1)):
        tmin = t
        t += 1
    return [[[bubble[col][t][x] for x in range(xmin, xmax+1)] for t in range(tmin, T)] for col in range(len(bubble))]

def center_bubble_again(bubble, kind):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    ext = int(N//8)
    if ext < 100: ext = 200
    if kind < 0:
        x1 = 0
        peaks1 = [phi_initial]
        while np.cos(np.nanmean(peaks1)) < np.cos(3.25) and x1 < N-1:
            slice = [bubble0[i][x1] for i in range(T)]
            peaks1 = [i for i in slice if i < mask]
            peaks2 = [i for i in range(T) if slice[i] < mask]
            x1 += 1
        if len(peaks2) > 0:    
            tmax = peaks2[-1]
        else:
            return 'next'
#        print('tmax, x1', tmax, x1)
        xmin = x1 - ext
        if xmin < 0: xmin = 0
        xmax = next((i for i in range(N-1, xmin, -1) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), N-1)
    else:
        x1 = N-1
        peaks1 = [phi_initial]
        while np.cos(np.nanmean(peaks1)) < np.cos(3.25) and x1 > 0:
            slice = [bubble0[i][x1] for i in range(T)]
            peaks1 = [i for i in slice if i < mask]
            peaks2 = [i for i in range(T) if slice[i] < mask]
            x1 -= 1
        if len(peaks2) > 0:    
            tmax = peaks2[-1]
        else:
            return 'next'
#        print('tmax, x1', tmax, x1)
        xmax = x1 + ext
        if xmax > N-1: xmax = N-1
        xmin = next((i for i in range(xmax) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), 0)

    t, tmin = 0, 0
    while not any([mask > bubble0[t][i] > phi_initial for i in range(xmin, xmax+1)]):
        tmin = t
        t += 1

    return [[[bubble[col][t][x] if bubble[col][t][x] != mask else normal[col] for x in range(xmin, xmax+1)] for t in range(tmin, tmax+1)] for col in range(len(bubble))]

def boost(vCOM, ga, x, t):
    return ga * (x - vCOM*t)

def interpolate_bubble(bubble, res):
    NT, N = len(bubble), len(bubble[0])
    bubble = np.asarray(bubble)
    t, x = np.arange(NT), np.arange(N)
    tt, xx = np.meshgrid(t, x)
    f = intp.interp2d(t, x, bubble[tt, xx], kind = 'quintic')
    T, X = np.arange(0, NT, 1/res), np.arange(0, N, 1/res)
    return f(T, X).T

def boost_bubble(bubble, vCOM, res):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    limit = 2*phi_initial
    ga = gamma(vCOM)

    t0 = int(time_at_size(bubble0, 100, limit)) # arbitrary choice of bubble centre
    x0 = np.mean([i for i in range(N) if mask > bubble0[t0][i] > limit])
    print('boost center: ', t0, x0)
    T0, X0 = boost(vCOM, ga, t0, x0), boost(vCOM, ga, x0, t0)
    deltaT, deltaX = np.abs(T0-t0), np.abs(X0-x0)

    bubble = [interpolate_bubble(bubble[col], res) for col in range(len(bubble))]
    new_bubble = [[[mask for x in range(N)] for t in range(T)] for col in range(len(bubble))]
    for t in range(T):
        for x in range(N):
            tt = int((boost(vCOM, ga, t, x) + np.sign(T//2-t0) * deltaT)*res)
            xx = int((boost(vCOM, ga, x, t) + np.sign(N//2-x0) * deltaX)*res)
            if (T*res > tt >= 0 and N*res > xx >= 0):
                for col in range(len(bubble)):
                    new_bubble[col][t][x] = bubble[col][tt][xx]
    return new_bubble, x0 - N//2

def the_whole_shebang(resolution, saveF):
    for sim in range(len(all_data))[::]:
#        exists = os.path.exists(bubbles_at_rest_file(sims_to_keep[sim]))
#        if not exists:
        if True:
            print(sim)
            bubble, dir = center_bubble(all_data[sim])
            beta = velocity(bubble[0], True)
            if beta == 'next':
                print('sim '+str(sim)+' skipped')
                continue
            elif np.abs(beta) > 0.9:
                beta = np.sign(beta)*0.9
            bubble = multiply_bubble(bubble, dir, fold(beta))

            fig = plt.figure(figsize=(5,5))
            plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
            plt.show()

            vtotal, bool1, bool2 = 0, True, False
            while np.abs(beta) > 0.5:
                bool2 = True
                bubble, dir = boost_bubble(bubble, beta, resolution)
                bubble = center_bubble_again(bubble, dir)
                if bubble == 'next':
                    print('sim '+str(sim)+' interrupted')
                    bool1 = False
                    break

                fig = plt.figure(figsize=(5,5))
                plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
                plt.show()

                prebeta = beta
                vtotal = add_velocities(vtotal, beta)
                print('step done ', beta)

                min_bubble = cut_out_bubble(bubble)
                fig = plt.figure(figsize=(5,5))
                plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
                plt.show()

                beta = velocity(min_bubble[0], True)
                if beta == 'next':
                    print('sim '+str(sim)+' interrupted')
                    if np.abs(prebeta) >= 0.9:
                        bool1 = False
                    break
                elif np.abs(beta) > 0.9:
                    beta = np.sign(beta)*0.9

            if not bool2:
                if np.abs(beta) > 0:
                    min_bubble = cut_out_bubble(bubble)
                    vtotal = beta
                    bool2 = True

            if bool1 and bool2:
                print('sim', sim, 'total vel ', vtotal)
                if saveF:
                    np.save(bubbles_at_rest_file(sims_to_keep[sim]), [min_bubble, vtotal])
    return

the_whole_shebang(3, True)
#the_whole_shebang(1, False)

0
1385
tcheckm, xmax, incr, decr  914 2048 20 29
tmin, tmax 644.0 2061.0
hyperbolic trajectories fit:  [2510.76244852  425.34454926 3337.20674657] [2108.03103795  838.83483207  730.40251877]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  986 2097.89898989899

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:30: RuntimeWarning: Mean of empty slice

sim 0 interrupted
1
-9
tcheckm, xmax, incr, decr  3064 2048 29 20
tmin, tmax 2834.0 4216.0
hyperbolic trajectories fit:  [1998.59245141 2994.79247779 2050.56979937] [1833.85872686 2824.81980186  197.03473364]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt
/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3098 6112.12

step done  -0.8609794053487859

tmin, tmax 1990.0 2431.0
hyperbolic trajectories fit:  [2181.82435298 2071.01062083 1078.12041181] [2094.75386882 1974.74215597 2137.98370655]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  2810 2176.8350515463917

step done  -0.7799861282372701

tmin, tmax 719.0 920.0

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

sim 1 interrupted
sim 1 total vel  -0.9817017601519805
2
1922
tcheckm, xmax, incr, decr  2920 2048 34 15
tmin, tmax 2761.0 4073.0
hyperbolic trajectories fit:  [1932.99051112 2718.50286694  407.13161464] [2059.72716534 2845.14279636  970.32472149]

boost center:  3013 6153.540816326531

step done  -0.8537547813578875

tmin, tmax 1933.0 2458.0
hyperbolic trajectories fit:  [2176.6265657  2082.72895618 1239.87604755] [ 1868.62549503  1780.94523772 10616.17003777]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  2723 2179.765306122449

step done  -0.9

tmin, tmax 406.0 564.0
sim 2 interrupted
3

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

-1795
tcheckm, xmax, incr, decr  3053 2048 28 21
tmin, tmax 2830.0 4222.0
hyperbolic trajectories fit:  [2092.35210468 2966.59332083 2891.98575649] [ 1884.88921763  2822.22167084 15943.46195861]

boost center:  3118 6163.13829787234

step done  -0.6971691359446275

tmin, tmax 2692.0 3410.0
hyperbolic trajectories fit:  [3848.08473389 2814.57405292 5332.16393492] [ 4065.41202846  2684.57837395 18896.75006022]
boost center:  3085 3887.7676767676767

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:46: RuntimeWarning: Mean of empty slice

sim 3 interrupted
4
96
tcheckm, xmax, incr, decr  3052 2048 20 29
tmin, tmax 2889.0 4214.0
hyperbolic trajectories fit:  [2149.61851055 2873.02520083  405.51417976] [2036.75153144 2987.24650157 1909.74219158]

boost center:  3137 2042.1546391752577

step done  0.8396396738711657

sim 4 interrupted
sim 4 total vel  0.8396396738711657
5
-1285
tcheckm, xmax, incr, decr  2600 2048 31 18
tmin, tmax 2460.0 3758.0
hyperbolic trajectories fit:  [2029.16230863 2527.25560879  573.67669492] [1920.1188301  2424.19013649 5500.41345676]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  2704 6099.010101010101

step done  -0.8043659632411795

tmin, tmax 2107.0 2614.0
hyperbolic trajectories fit:  [2395.26959571 2166.14145269 2839.45758867] [2480.82044105 2113.19387923 4105.91900267]
sim 5 total vel  -0.8043659632411795
6
544
tcheckm, xmax, incr, decr  3068 2048 38 11
tmin, tmax 2733.0 4223.0
hyperbolic trajectories fit:  [1717.10893247 2625.83785454  322.42848807] [2083.20683125 2998.54863118 1073.72224572]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  3127 6169.020202020202

step done  -0.9

tmin, tmax 1673.0 2114.0
hyperbolic trajectories fit:  [1392.615418   1802.91944168 1814.69911372] [1167.51123732 1575.79395864 1702.16869512]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  2599 1792.157894736842

step done  -0.9

tmin, tmax 238.0 403.0
sim 6 interrupted
7
-1194
tcheckm, xmax, incr, decr  3047 2048 22 27
tmin, tmax 2831.0 4202.0
hyperbolic trajectories fit:  [2043.78271308 2972.40988835 3437.04932082] [2208.37372309 2816.64471676 3671.56005041]

boost center:  3150 2060.87

step done  0.7279749996226531

sim 7 interrupted
sim 7 total vel  0.7279749996226531
8
-1981
tcheckm, xmax, incr, decr  2097 2048 15 34
tmin, tmax 1927.0 3260.0
hyperbolic trajectories fit:  [2085.24669814 2031.94404518 1911.4220065 ] [2229.30213412 1888.41724278 2734.12620153]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  2185 2099.776595744681

step done  0.8165919266899307

tmin, tmax 845.0 847.0
sim 8 interrupted
sim 8 total vel  0.8165919266899307

/home/dpirvu/.local/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/home/dpirvu/.local/lib/python3.7/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

9
965
tcheckm, xmax, incr, decr  3056 2048 14 35
tmin, tmax 2861.0 4217.0
hyperbolic trajectories fit:  [2024.24737224 2988.34513192 2183.47048059] [2254.76148444 2760.18547405  334.47970235]

boost center:  3143 2059.81

step done  0.9

tmin, tmax 1697.0 2142.0
hyperbolic trajectories fit:  [ 637.90854895 1831.13015548  634.79033482] [1003.835134   1460.87103372 7795.97776206]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  2524 1269.07

step done  0.9

tmin, tmax 339.0 467.0
sim 9 interrupted
10
-571
tcheckm, xmax, incr, decr  3063 2048 29 20
tmin, tmax 2867.0 4218.0
hyperbolic trajectories fit:  [1933.21636791 2855.52827496  515.97330125] [2070.53348553 2988.07900252 1960.9483611 ]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3072 6156.62

step done  -0.8433814632785198

tmin, tmax 2118.0 2570.0
hyperbolic trajectories fit:  [2323.69972345 2128.56518008 -206.01188658] [2352.45902813 2189.9203633   297.89443905]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt
/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt
/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:78: RuntimeWarning: invalid value encountered in sqrt
/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:79: RuntimeWarning: invalid value encountered in sqrt

boost center:  2890 2379.3125

step done  -0.7845853874876656

tmin, tmax 737.0 934.0
sim 10 interrupted
sim 10 total vel  -0.9796968017633337

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

11
819
tcheckm, xmax, incr, decr  2246 2048 18 31
tmin, tmax 2111.0 3406.0
hyperbolic trajectories fit:  [2139.62035895 2073.78038008  201.13576774] [2037.3809434  2180.5088944  1237.49791624]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  2337 2014.7448979591836

step done  0.8748090434334704

sim 11 interrupted
sim 11 total vel  0.8748090434334704
12
718
tcheckm, xmax, incr, decr  3057 2048 25 24
tmin, tmax 2864.0 4222.0
hyperbolic trajectories fit:  [2067.90664455 2983.6513148  2000.08795986] [1888.89367526 2833.17572067 9980.51451464]

boost center:  3144 6163.79797979798

step done  -0.764586479591757

tmin, tmax 2541.0 3142.0
hyperbolic trajectories fit:  [ 3311.75462531  2518.78367497 13702.76225382] [3141.45024919 2637.84924056 2800.01462095]
sim 12 total vel  -0.764586479591757
13
201
tcheckm, xmax, incr, decr  3139 2048 24 25
tmin, tmax 2875.0 4217.0
hyperbolic trajectories fit:  [1929.5775653  2990.04209954 1358.09809565] [1719.6253256  2780.48066199 5777.30554417]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  3149 1949.58

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:30: RuntimeWarning: Mean of empty slice

sim 13 interrupted
14
-561
tcheckm, xmax, incr, decr  3056 2048 21 28
tmin, tmax 2932.0 4195.0
hyperbolic trajectories fit:  [2008.91105665 2917.99472962 1553.45167839] [2051.92842748 2968.63521851   21.16504674]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

sim 14 total vel  -0.46042050198857987
15
1967
tcheckm, xmax, incr, decr  1783 2048 17 32
tmin, tmax 1594.0 2943.0
hyperbolic trajectories fit:  [2077.11990151 1713.0985189  2120.38145321] [2296.45770193 1499.68014336 3550.78483524]

boost center:  1856 2068.021276595745

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:30: RuntimeWarning: Mean of empty slice

sim 15 interrupted
16
1797
tcheckm, xmax, incr, decr  3060 2048 40 9
tmin, tmax 2910.0 4222.0
hyperbolic trajectories fit:  [2050.59960987 2995.27354738 1594.98593567] [1835.03241416 2777.9706965  3877.36549426]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

sim 16 total vel  0.13951297613441357
17
-852
tcheckm, xmax, incr, decr  3062 2048 29 20
tmin, tmax 2948.0 4220.0
hyperbolic trajectories fit:  [1879.9872072  2866.48690605 4865.44905482] [2002.4877874  2998.04285695 1196.65731133]

boost center:  3131 6099.944444444444

step done  -0.6019582065205895

tmin, tmax 2891.0 3610.0
hyperbolic trajectories fit:  [5077.40902223 2932.53981519  477.0458381 ] [ 4938.5509182   2857.76325842 17040.83380892]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  3033 5058.727272727273

step done  -0.6779969034076327

tmin, tmax 1636.0 2046.0
hyperbolic trajectories fit:  [ 1575.80792508  1657.83907576 35567.60548812] [ 2069.77073255  1388.23117192 33976.18665489]
boost center:  3048 1728.3367346938776

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:46: RuntimeWarning: Mean of empty slice

sim 17 interrupted
18
569
tcheckm, xmax, incr, decr  3062 2048 22 27
tmin, tmax 2788.0 4216.0
hyperbolic trajectories fit:  [2256.90579656 2790.89661684  183.78571978] [2065.59119416 2995.70344859 1069.49406554]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3070 2056.94

step done  0.9

tmin, tmax 1515.0 1665.0
sim 18 interrupted
19
-1518
tcheckm, xmax, incr, decr  3063 2048 22 27
tmin, tmax 2884.0 4225.0
hyperbolic trajectories fit:  [2028.25058499 2981.85036645 3018.36625944] [2142.65446657 2901.15849357 4802.28853571]

boost center:  3138 2061.5858585858587

step done  0.5286828593058204

sim 19 interrupted
sim 19 total vel  0.5286828593058204
20
1302
tcheckm, xmax, incr, decr  3057 2048 28 21
tmin, tmax 2813.0 4215.0
hyperbolic trajectories fit:  [2057.45266121 2988.19270195 2003.98206648] [ -2040.31605517  -1114.32063744 -36713.88378622]

sim 20 total vel  -0.006909779998893939
21
148
tcheckm, xmax, incr, decr  3069 2048 33 16
tmin, tmax 2770.0 4221.0
hyperbolic trajectories fit:  [1742.78465684 2696.51672289  624.20108106] [2038.77104365 2997.59020138 1454.75184529]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  3150 6144.15

step done  -0.9

tmin, tmax 1691.0 2106.0
hyperbolic trajectories fit:  [1638.80395057 1647.03891516 1522.88547097] [1785.32453116 1800.63913954 1148.16631298]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  2533 1787.0

step done  -0.7760460228056966

tmin, tmax 591.0 771.0
sim 21 interrupted
sim 21 total vel  -0.9868141476951815
22
1565
tcheckm, xmax, incr, decr  1090 2048 22 27
tmin, tmax 945.0 2246.0
hyperbolic trajectories fit:  [2011.85020529 1018.32650523 1509.66052844] [2124.45881166  906.83071679 2295.55550735]

boost center:  1097 2021.57

step done  0.7860618056377217

tmin, tmax 908.0 1403.0
hyperbolic trajectories fit:  [ 706.61383561  919.5571983  4001.11409413] [ 779.49734586  929.09147054 1473.91163883]
sim 22 total vel  0.7860618056377217
23
-887
tcheckm, xmax, incr, decr  1673 2048 26 23
tmin, tmax 1453.0 2762.0
hyperbolic trajectories fit:  [ 1715.01809593  1095.31540593 12437.15336213] [2157.5446663  1534.6407919  1368.01613859]

boost center:  1694 6275.33

step done  -0.7761008773387426

tmin, tmax 1409.0 1918.0
hyperbolic trajectories fit:  [  714.38097504  1320.27972241 14241.52313717] [ 866.01994562 1458.52294677  793.57661621]
sim 23 total vel  -0.7761008773387426

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

24
907
tcheckm, xmax, incr, decr  2178 2048 17 32
tmin, tmax 2023.0 3334.0
hyperbolic trajectories fit:  [2190.14919044 1976.42349681  645.51398929] [2060.86887774 2105.01153307 1201.50201722]

boost center:  2186 2061.0714285714284

step done  0.8217883904067935

sim 24 interrupted
sim 24 total vel  0.8217883904067935
25
1478
tcheckm, xmax, incr, decr  3059 2048 17 32
tmin, tmax 2855.0 3944.0
hyperbolic trajectories fit:  [2374.34422393 2672.77012558  996.22936561] [2060.48465552 2991.5316991   869.06708118]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3069 2045.080808080808

step done  0.9

tmin, tmax 1860.0 1985.0
sim 25 interrupted
26

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

1074
tcheckm, xmax, incr, decr  3069 2048 25 24
tmin, tmax 2928.0 4222.0
hyperbolic trajectories fit:  [2045.95093166 2984.38196391 1965.81494544] [1969.35622525 2928.44808831 3521.88604484]

sim 26 total vel  -0.4828900790231301
27
-23
tcheckm, xmax, incr, decr  3062 2048 31 18
tmin, tmax 2918.0 4211.0
hyperbolic trajectories fit:  [2015.11530867 2990.48936496 1042.49806198] [1948.75356718 2912.33148266  242.59700974]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

boost center:  3066 6120.989795918367

step done  -0.740797160006187

tmin, tmax 2691.0 3258.0
hyperbolic trajectories fit:  [3422.40512129 2725.43445155 3838.40500812] [3357.50802008 2707.78944152 1043.44051934]
sim 27 total vel  -0.740797160006187

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in sqrt

28
1113
tcheckm, xmax, incr, decr  3068 2048 19 30
tmin, tmax 2828.0 4223.0
hyperbolic trajectories fit:  [2385.95127705 2656.33815381 1982.85265188] [2051.53819468 2998.12552642 1079.88073432]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3080 2035.4421052631578

step done  0.8785624964395339

sim 28 interrupted
sim 28 total vel  0.8785624964395339
29
-224
tcheckm, xmax, incr, decr  3066 2048 32 17
tmin, tmax 2882.0 4228.0
hyperbolic trajectories fit:  [2062.34367886 2986.72521452 3302.72545845] [1925.14103226 2883.1849998  5449.1188915 ]

boost center:  3089 6143.818181818182

step done  -0.5883580653525555

tmin, tmax 2849.0 3650.0
hyperbolic trajectories fit:  [5308.81903084 2937.10004699 4931.607496  ] [ 5103.56502234  2852.04640042 22887.56255026]
sim 29 total vel  -0.5883580653525555
30
316
tcheckm, xmax, incr, decr  3060 2048 16 33
tmin, tmax 2919.0 4213.0
hyperbolic trajectories fit:  [2142.32151983 2878.02662749  251.60951299] [2037.69467296 2988.62320562 1006.5062018 ]

boost center:  3127 2019.141414141414

step done  0.8657470271084928

sim 30 interrupted
sim 30 total vel  0.8657470271084928
31
1504
tcheckm, xmax, incr, decr  3067 2048 17 32
tmin, tmax 2730.0 4219.0
hyperbolic trajectories fit:  [2402.530071   2636.83761415  135.89361914] [2049.98950266 2993.22408142 1184.60547628]

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

boost center:  3075 2041.4141414141413

step done  0.9

tmin, tmax 1544.0 1963.0
hyperbolic trajectories fit:  [ 616.04011905 1722.29546096 1676.56720445] [  2974.82564624   -612.35750717 124226.05845007]
sim 31 total vel  0.9

/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in sqrt

32
1182
tcheckm, xmax, incr, decr  3147 2048 29 20
tmin, tmax 2898.0 4218.0
hyperbolic trajectories fit:  [2239.74398221 2941.94163778 4187.65406162] [2130.91996524 2934.05271353 5850.68557122]

sim 32 total vel  -0.060788704260273003
33
336
tcheckm, xmax, incr, decr  3071 2048 20 29
tmin, tmax 2852.0 4232.0
hyperbolic trajectories fit:  [2029.23861921 2981.22811141 3495.26486516] [ 2213.72792097  2851.96513508 11736.5028761 ]

boost center:  3159 2074.8125

step done  0.6422838051923997

sim 33 interrupted
sim 33 total vel  0.6422838051923997

Average bubble

 def cut_out_coordinates(bubble0):
    T, N = len(bubble0), len(bubble0[0])
    limit = np.ceil(phi_upper_bound)
    size = 120 # do not change

    tmin = T-1
    for t in range(T):
        slice = smoothen(bubble0[t], dx, False)
        right_phi = [x for x in range(N-1) if (slice[x] >= limit and slice[x+1] >= limit)]
        if len(right_phi) >= size:
            tmin = t
            break
    if tmin == T-1:
        return None
    xcen = int(np.round(np.nanmean(right_phi)))
    xcen_list = [np.round(np.nanmean([x for x in range(N) if bubble0[t][x] >= limit])) for t in range(tmin-50,tmin,2)]
    if stat.stdev(xcen_list) > 50:
        return None

#    fig = plt.figure(figsize=(4,4))
#    plt.plot(xcen, tmin, 'ro')
#    plt.imshow(bubble, aspect='auto', interpolation='none', origin='lower')
#    plt.show()
    ex1, ex2, ex3, ex4 = np.arange(tmin)[::-1], np.arange(tmin, T), np.arange(xcen)[::-1], np.arange(xcen, N)

    if len(ex1) > 500: ex1 = ex1[:500]
    if len(ex2) > 500: ex2 = ex2[:500]
    if len(ex3) > 500: ex3 = ex3[:500]
    if len(ex4) > 500: ex4 = ex4[:500]    
    return ex1, ex2, ex3, ex4

def stack_bubbles(bubble_data):
    extent_list, bad_bubbles = [], []

    for sim in range(len(bubble_data)):
        coords = cut_out_coordinates(bubble_data[sim][0])
        if coords is not None:
            extent_list.append(coords)
        else:
            bad_bubbles.append(sim)

    extent_list = np.asarray(extent_list)
    print('bad bubbles', bad_bubbles)

    bubble_data = [bubble_data[i] for i in range(len(bubble_data)) if i not in bad_bubbles]
    nCols = 5
    for sim in range(len(bubble_data)):
        bubble_data[sim].append([0.5*np.asarray(i)**2 for i in bubble_data[sim][1]])

    uper_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    upper_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    lower_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    lower_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]

    av_mat, av_err_mat = [], []
    for col in range(nCols):
        average_matrix, average_error_matrix = [[],[],[],[]], [[],[],[],[]]
        coord_lists = [uper_right_bubble_list[col], upper_left_bubble_list[col], lower_right_bubble_list[col], lower_left_bubble_list[col]]

        for i in range(4):
            max_lines, max_cols = 500, 500
            for line in range(max_lines):
                average_matrix[i].append([normal[col] for n in range(max_cols)])
                average_error_matrix[i].append([0. for n in range(max_cols)])

            for num_line in range(max_lines):
                for num_col in range(max_cols):
                    meas = []
                    for simulation in coord_lists[i]:
                        if len(simulation) > num_line and len(simulation[0]) > num_col:
                            meas.append(simulation[num_line][num_col])

                    average_matrix[i][num_line][num_col] = np.mean(meas)
                    if len(meas) > 1.:
                        average_error_matrix[i][num_line][num_col] = stat.stdev(meas)/len(meas)
        
        av_mat.append(average_matrix)
        av_err_mat.append(average_error_matrix)
    return [av_mat, av_err_mat], bad_bubbles

def restore_average_bubble(stacked_bubble_data):
    whole_bubbles = []
    for av_bub in stacked_bubble_data:
        whole_bubbles.append([])
        for col in range(len(av_bub)):
            top = np.concatenate((np.flip(np.asarray(np.flip(av_bub[col][0],1)).transpose(),1), np.flip(np.asarray(av_bub[col][1]).transpose(),0)), axis=1)
            bottom = np.concatenate((np.flip(np.asarray(av_bub[col][2]).transpose(),1), np.flip(np.asarray(np.flip(av_bub[col][3],0)).transpose(),1)), axis=1)
            whole_bubbles[-1].append(np.concatenate((top, bottom), axis=0))
    return np.asarray(whole_bubbles)

bubble_list = []
for sim in range(nSims):
    bool = os.path.exists(bubbles_at_rest_file(sim))
    if bool:
        bubble_list.append(np.load(bubbles_at_rest_file(sim))) #deboosted_bubbles2 #deboosted_bubbles_loop_all_sims7
bubble_list = np.asarray(bubble_list)
bubble_list_save = bubble_list[:,0]
bubble_velocities = bubble_list[:,1]

#if testing
#normal = [phi_initial, 0, 0, V(phi_initial), 0]

#either compute average
stacked_bubble_data, bad_bubbles = stack_bubbles(bubble_list_save)
ab, error_ab = restore_average_bubble(stacked_bubble_data)
np.save(average_bubble(minSim, nSims), [ab, error_ab])
#or load it
#ab, error_ab = np.load(average_bubble(minSim, nSims))

#separate columns
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
error_bubble, error_mom_bubble, error_ge_bubble, error_pe_bubble, error_ke_bubble = error_ab

bad bubbles [1, 3, 4, 6, 10, 11, 16, 19, 21, 24]

#plot discarded bubbles
for i in range(len(bubble_list_save)):
    if i in bad_bubbles:
        fig, ax = plt.subplots(1, 1, figsize=(7,5))
        im = ax.imshow(bubble_list_save[i][0], aspect='auto', interpolation='none', origin='lower')
        clb = plt.colorbar(im, ax = ax)
        plt.show()

bubble_list_save = [bubble_list_save[i] for i in range(len(bubble_list_save)) if i not in bad_bubbles]
bubble_velocities = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i not in bad_bubbles]
check = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i in bad_bubbles]
print('Total number of bubbles available: ', len(bubble_list_save))

Total number of bubbles available:  15

for i in bubble_list_save:
    fig, ax = plt.subplots(1, len(i), figsize=(25,5))
    for col in range(len(i)):
        im = ax[col].imshow(i[col], aspect='auto', interpolation='none', origin='lower')
        clb = plt.colorbar(im, ax = ax[col])
    plt.show()

#plot average bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble
fig, ax = plt.subplots(1, len(ab), figsize=(25, 5))
for col in range(len(ab)):
    im = ax[col].imshow(ab[col], aspect='auto', interpolation='none', origin='lower')
    clb = plt.colorbar(im, ax = ax[col])
plt.savefig(plots_file+'average_bubble'+suffix+'.png')

#plot distribution of boosts
fig, ax = plt.subplots(1, 3, figsize=(15,5))
ax[0].hist(bubble_velocities, bins=len(bubble_velocities))
ax[1].hist([gamma(i) for i in bubble_velocities], bins=len(bubble_velocities))
ax[2].hist([rapidity(i) for i in bubble_velocities], bins=len(bubble_velocities))

[ax[i].set_ylabel('Count') for i in range(len(ax))]
ax[0].set_xlabel(r'$v$'); ax[1].set_xlabel(r'$\gamma$'); ax[2].set_xlabel(r'$w$');
plt.savefig(plots_file+'boost_distrib'+suffix+'.png');

def coleman_match(bubble, colemanSoln):
    minSum = np.sum(colemanSoln)**2.
    tColeman = 0
    for t in range(len(bubble)):
        slice = smoothen(bubble[t], 2*dx, True)
        a = sum([(slice[x] - colemanSoln[x])**2. for x in range(len(slice))])
        if a < minSum:
            minSum, tColeman = a, t
    return tColeman

def plot_coleman_match(bubble, error_bubble, colemanSoln):
    T, N = len(bubble), len(bubble[0])
    colemanSoln = colemanSoln[len(colemanSoln)//2-N//2:len(colemanSoln)//2+N//2]
    tcoleman = coleman_match(bubble, colemanSoln)
    slice, err_slice = bubble[tcoleman], error_bubble[tcoleman]

    fig = plt.figure(figsize=(15,5))
    plt.errorbar(np.arange(N), slice, yerr = err_slice, color = 'k', marker='o', ms=1, alpha = 0.8, ecolor = 'green', label = 'Average bubble')
    plt.plot(np.arange(N), colemanSoln, 'r-', label = 'Best match t = '+str(tcoleman)); plt.legend()
    plt.savefig(plots_file+'average_bubble_profile'+suffix+'.png')
    return tcoleman

def bubble_total_energy(KED, GED, PED):
    T, N = len(KED), len(KED[0])
    KE = [sum(KED[t]) for t in range(T)]
    GE = [sum(GED[t]) for t in range(T)]
    PE = [sum(PED[t]) for t in range(T)]
    TE = [sum(x) for x in zip(KE, GE, PE)]
    return KE, GE, PE, TE

def plot_total_energy_density(KED, GED, PED, tcoleman):
    energies = [KED, GED, PED, KED + GED + PED]
    titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']

    fig, ax = plt.subplots(2, 2, figsize=(15,10))
    N = len(KED[0])
    for i in range(4):
        if i in [0,1]: j = 0
        else: j = 1
        im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower') 
        clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
        ax[i%2][j].plot(np.arange(N), np.ones(N)*tcoleman, 'r-')
    plt.savefig(plots_file+'energy_densities_from_average'+suffix+'.png');
    return

def plot_total_energy(KED, GED, PED, tcoleman):
    energies = bubble_total_energy(KED, GED, PED)
    titles = ['KE', 'GE', 'PE', 'TE']
    ls = ['b-', 'g-', 'y-', 'r-']
    N = len(energies[0])
    fig = plt.figure(figsize=(8,5))
    for i in range(len(energies)):
        plt.plot(np.arange(N), energies[i], ls[i], label = titles[i])
    plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
    plt.axvline(tcoleman, color = 'darkgray', ls = '--')
    plt.xlabel('$t$'); plt.legend(); plt.savefig(plots_file+'energies'+suffix+'.png'); plt.show()
    return

def plot_energy_dens_at_t(KED, GED, PED, timeslice):
    TED = [sum(x) for x in zip(KED, GED, PED)]
    energies = [KED, GED, PED, TED]
    titles = ['KED', 'GED', 'PED', 'TED']
    ls = ['b-', 'g-', 'y-', 'r-']
    N = len(KED[0])
    
    fig = plt.figure(figsize=(8,8)); plt.subplots_adjust(hspace=.0)
    gs = gridspec.GridSpec(len(energies), 1, height_ratios=[1]*len(energies))
    ax = [[]]*len(energies); ax[0] = plt.subplot(gs[0])
    for i in range(1, len(ax)): ax[i] = plt.subplot(gs[i], sharex = ax[i-1])
    for i in range(len(ax)):
        ax[i].plot(np.arange(N), energies[i][timeslice], ls[i], label = titles[i])
        if i == 0:
            ax[i].axhline(np.pi, color = 'darkgray', ls='--')
        else:
            ax[i].axhline(np.mean(energies[i][timeslice]), color = 'darkgray', ls='--')            
        ax[i].legend()
    ax[-1].set(xlabel='$x$')
    plt.savefig(plots_file+'energy_densities_at_t'+str(timeslice)+suffix+'.png'); plt.show()
    return

TColeman = plot_coleman_match(bubble, error_bubble, coleman_profile)
plot_total_energy_density(ke_bubble, ge_bubble, pe_bubble, TColeman)
plot_total_energy(ke_bubble, ge_bubble, pe_bubble, TColeman)
#plot_energy_dens_at_t(ke_bubble, ge_bubble, pe_bubble, TColeman)

def GradEnDen(field):
    Nx = len(field)
    NNx = Nx//2+1
    fft_field = scp.fft.fft(field)
    for i in range(len(fft_field)):
        fft_field[i] = 1j*i*dk * fft_field[i]# / Nx
    real_field = scp.fft.ifft(fft_field)
    return [0.5*i.real**2 for i in real_field]

def plot_derived_total_energy_density(field, momentum, filter, t_coleman):
    smooth_field = [smoothen(slice, filter, True) for slice in field]
    smooth_momentum = [smoothen(slice, filter, True) for slice in momentum]
    KED = [0.5*slice**2. for slice in smooth_momentum]
    GED = [GradEnDen(slice) for slice in smooth_field]
    PED = [V(slice) for slice in smooth_field]
    TED = [[KED[t][x]+GED[t][x]+PED[t][x] for x in range(len(KED[0]))] for t in range(len(KED))]
    energies = [KED, GED, PED, TED]
    titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']

    fig, ax = plt.subplots(2, 2, figsize=(15,10))
    N = len(smooth_field[0])
    for i in range(4):
        if i in [0,1]: j = 0
        else: j = 1
        im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower') 
        clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
        ax[i%2][j].plot(np.arange(N), np.ones(N)*t_coleman, 'r-')
    plt.savefig(plots_file+'energy_densities'+suffix+'.png');
    return

def plot_TE_with_Nbubbles(tcoleman):
    bubble_averages_list = []
    for i in range(Nbubs):
        bool = os.path.exists(average_of_N_bubbles(i+1))
        if bool:
            bubble_averages_list.append(np.load(average_of_N_bubbles(i+1)))

    fig = plt.figure(figsize=(8,6))
    for bubble_average in bubble_averages_list:
        ab, error_ab, length = bubble_average
        bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
        energies = bubble_total_energy(ke_bubble, ge_bubble, pe_bubble)
        plt.plot(np.arange(len(energies[-1])), energies[-1], label = length)
        plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
    plt.axvline(tcoleman, color = 'darkgray', ls = '--')
    plt.xlabel('$t$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_nbubbles1'+suffix+'.png');

    fig = plt.figure(figsize=(8,6))
    for bubble_average in bubble_averages_list:
        ab, error_ab, length = bubble_average
        bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
        GED, PED, KED = ge_bubble[tcoleman], pe_bubble[tcoleman], ke_bubble[tcoleman]
        TED = GED + PED + KED
        plt.plot(length, sum(TED), 'ro')#, np.mean(TED))
    plt.xlabel('$N$'); plt.ylabel('$E$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_nbubbles2'+suffix+'.png');
    return

def plot_TE_with_filter(tcoleman):
    bubble_average = np.load(average_of_N_bubbles(Nbubs))
    ab, error_ab, length = bubble_average
    bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
    slice = bubble[tcoleman]
    mom_slice = mom_bubble[tcoleman]

    fig = plt.figure(figsize=(8,6))
    for filter in np.linspace(dx, 1000*dx, 200):
        smooth_field = smoothen(slice, filter, True)
        smooth_momentum = smoothen(mom_slice, filter, True)
        KED = 0.5*smooth_momentum**2.
        GED = GradEnDen(smooth_field)
        PED = V(smooth_field)
        TED = KED + GED + PED
        plt.plot(filter, np.sum(TED), 'ro', ms=2)#, 'o', label='{:.3f}'.format(filter))
    plt.axhline(0, ls='--', color='darkgray')
    plt.xlabel(r'$\sigma$'); plt.ylabel('$E$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_filter'+suffix+'.png');
    return

plot_derived_total_energy_density(bubble, mom_bubble, 1, TColeman)

#average x bubbles at a time; save
Nbubs = len(bubble_list_save)
step = 2
for i in range(0, Nbubs, step):
    qq = len(bubble_list_save)
    bubble_puzzle, _ = stack_bubbles(bubble_list_save)
    ab, error_ab = restore_average_bubble(bubble_puzzle)
    np.save(average_of_N_bubbles(qq), [ab, error_ab, qq])
    for i in range(step):
        if len(bubble_list_save) > 1:
            del bubble_list_save[0]
        else:
            break

bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []

plot_TE_with_Nbubbles(TColeman)
plot_TE_with_filter(TColeman)

No handles with labels found to put in legend.
No handles with labels found to put in legend.